home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / poly / balance.c next >
Encoding:
C/C++ Source or Header  |  2000-05-05  |  2.6 KB  |  128 lines
  1. /* poly/balance.c
  2.  * 
  3.  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Brian Gough
  4.  * 
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 2 of the License, or (at
  8.  * your option) any later version.
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20. static void balance_companion_matrix (double *m, size_t n);
  21.  
  22. #define RADIX 2
  23. #define RADIX2 (RADIX*RADIX)
  24.  
  25. static void
  26. balance_companion_matrix (double *m, size_t nc)
  27. {
  28.   int not_converged = 1;
  29.  
  30.   double row_norm = 0;
  31.   double col_norm = 0;
  32.  
  33.   while (not_converged)
  34.     {
  35.       size_t i, j;
  36.       double g, f, s;
  37.  
  38.       not_converged = 0;
  39.  
  40.       for (i = 0; i < nc; i++)
  41.     {
  42.       /* column norm, excluding the diagonal */
  43.  
  44.       if (i != nc - 1)
  45.         {
  46.           col_norm = fabs (MAT (m, i + 1, i, nc));
  47.         }
  48.       else
  49.         {
  50.           col_norm = 0;
  51.  
  52.           for (j = 0; j < nc - 1; j++)
  53.         {
  54.           col_norm += fabs (MAT (m, j, nc - 1, nc));
  55.         }
  56.         }
  57.  
  58.       /* row norm, excluding the diagonal */
  59.  
  60.       if (i == 0)
  61.         {
  62.           row_norm = fabs (MAT (m, 0, nc - 1, nc));
  63.         }
  64.       else if (i == nc - 1)
  65.         {
  66.           row_norm = fabs (MAT (m, i, i - 1, nc));
  67.         }
  68.       else
  69.         {
  70.           row_norm = (fabs (MAT (m, i, i - 1, nc)) 
  71.                           + fabs (MAT (m, i, nc - 1, nc)));
  72.         }
  73.  
  74.       if (col_norm == 0 || row_norm == 0)
  75.         {
  76.           continue;
  77.         }
  78.  
  79.       g = row_norm / RADIX;
  80.       f = 1;
  81.       s = col_norm + row_norm;
  82.  
  83.       while (col_norm < g)
  84.         {
  85.           f *= RADIX;
  86.           col_norm *= RADIX2;
  87.         }
  88.  
  89.       g = row_norm * RADIX;
  90.  
  91.       while (col_norm > g)
  92.         {
  93.           f /= RADIX;
  94.           col_norm /= RADIX2;
  95.         }
  96.  
  97.       if ((row_norm + col_norm) < 0.95 * s * f)
  98.         {
  99.           not_converged = 1;
  100.  
  101.           g = 1 / f;
  102.  
  103.           if (i == 0)
  104.         {
  105.           MAT (m, 0, nc - 1, nc) *= g;
  106.         }
  107.           else
  108.         {
  109.           MAT (m, i, i - 1, nc) *= g;
  110.           MAT (m, i, nc - 1, nc) *= g;
  111.         }
  112.  
  113.           if (i == nc - 1)
  114.         {
  115.           for (j = 0; j < nc; j++)
  116.             {
  117.               MAT (m, j, i, nc) *= f;
  118.             }
  119.         }
  120.           else
  121.         {
  122.           MAT (m, i + 1, i, nc) *= f;
  123.         }
  124.         }
  125.     }
  126.     }
  127. }
  128.